use ${newdata}carownership_dataset_bergen, clear
/* Only keeping necessary variables */
drop if $trmgroup == . | $trmgroup == 0
capt drop bergen 

egen long famid_num = group(familienr)
gen commuters = ($trmgroup == 1 | $trmgroup == 3)
gen bergen = ($trmgroup == 1 | $trmgroup == 2)
gen bergen_commuters = bergen * commuters

/* All time-varying variables based on geography (residence or work location)
   are set to their 2014 values (variation post 2014 is endogenous to treatment) */
foreach var in dist time_work PublicVSCarTime_fam_mean PublicDiffCarTime_fam_mean grk {
	replace `var' = . if year != 2014
	bysort familienr (`var'): replace `var' = `var'[_n-1] if missing(`var')
}


local years
local years_commuters
local years_grk
local years_bergen
local years_bergen_commuters
local years_stavanger_commuters
forvalues i = $firstyearpre/$lastyearpost {
	if `i' != $lastyearpre qui gen year_`i' = (year == `i')
	if `i' != $lastyearpre local years `years' year_`i'
	if `i' != $lastyearpre qui gen commuters_`i' = commuters * year_`i'
	if `i' != $lastyearpre local years_commuters `years_commuters' commuters_`i'
	if `i' != $lastyearpre local years_grk `years_grk' c.year_`i'#i.grk
	if `i' != $lastyearpre qui gen bergen_commuters_`i' = bergen_commuters * year_`i'
	if `i' != $lastyearpre local years_bergen_commuters `years_bergen_commuters' bergen_commuters_`i'
}



/* === Only the triple difference ===================*/
reghdfe $yvar ///
	`years_bergen_commuters' /// DiDiD: post commuter in bergen
	`years_commuters' /// 
	$xvar /// 
	commuters /// 
	bergen_commuters /// 
	, absorb( ///
		i.grk ///
		`years_grk' ///
		i.famid_num  ///
	) vce(cluster i.grk)	
estimates save "${ster}fig2${yvar}", replace			

/* Saving estimates */
matrix B = e(b)
matrix V = e(V)

/* Storing them as variables */
qui gen beta = .
qui gen se = .
qui gen xaxis = .
local coef = $lastyearpost - $firstyearpre /* The number of coef we want to keep */
forvalues i = 1/`coef' {
	qui replace xaxis = $firstyearpre + `i' - 1 in `i'
	qui replace beta  = B[1,`i'] in `i'
	qui replace se    = sqrt(V[`i',`i']) in `i'
}
/* Adding the base year here */
qui replace xaxis = xaxis + 1 		if xaxis >= $lastyearpre /* To skip the base year */	
local lastobs = `coef' + 1 			/* Adding an observation for the base year */
qui replace beta = 0 				in `lastobs' /* Coef is normalized to zero */
qui replace se = 	0 				in `lastobs'
qui replace xaxis = $lastyearpre 	in `lastobs'
/* 95 % CI */
qui gen plusse = beta + 1.96 * se
qui gen minusse = beta - 1.96 * se
gen time = mdy(1,1,xaxis)
format time %td

/* === FIGURE ============================================================*/
gen ttime = year(time)
local announcement_date = 2014 + (mdy(2,18,2015) - mdy(1,1,2015) + 1) / 365
local implementation_date = 2015 + (mdy(2,1,2016) - mdy(1,1,2016) + 1) / 366
if "$yvar" == "bev" {
	local ylabel "-0.02(0.02)0.08"
}
if "$yvar" == "ice" {
	local ylabel "-0.08(0.02)0.08"
}
if "$yvar" == "cars" {
	local ylabel "-0.08(0.02)0.08"
}
twoway ///	
	(rarea plusse minusse ttime ///
		, sort color(gs10%50)) ///			
	(scatter beta ttime ///
		, sort connect(l) lcolor(black) mcolor(black) lpattern(solid)) ///	
	, graphregion(color(white))	///
	xtitle("") ytitle("Triple difference estimate") ///
	name(fig4${yvar}, replace) ///
	legend(off order(2 1) ///
		label(1 "95 % CI") ///
		label(2 "Estimate") ///
	) ///
	xline(`announcement_date', lpattern(shortdash) lcolor(black)) ///mdy(2,1,2016) - implementation date
	xline(`implementation_date', lpattern(dash) lcolor(black)) /// /mdy(2,18,2015) - announcement date
	yline(0) ylabel(`ylabel') ///
	xlabel(2011(1)2017)

graph export "${figures}fig4${yvar}.png", as(png) replace
graph export "${figures}fig4${yvar}.pdf", as(pdf) replace
	


if "$yvar" == "bev" {
	gen outcome_bergen_commuter = .
	forvalues tt = $firstyearpre(1)$lastyearpost {
		sum $yvar if e(sample) == 1 & bergen_commuters == 1 & year == `tt'
		replace outcome_bergen_commuter = r(mean) if xaxis == `tt'
	}
	gen outcome_bergen_commuter2 = outcome_bergen_commuter - beta
	replace outcome_bergen_commuter2 = . if xaxis <= $lastyearpre
	replace outcome_bergen_commuter2 = outcome_bergen_commuter if xaxis == $lastyearpre

	replace outcome_bergen_commuter = outcome_bergen_commuter * 100	// measure in percent
	replace outcome_bergen_commuter2 = outcome_bergen_commuter2 * 100 // measure in percent

	
	local announcement_date = 2014 + (mdy(2,18,2015) - mdy(1,1,2015) + 1) / 365
	local implementation_date = 2015 + (mdy(2,1,2016) - mdy(1,1,2016) + 1) / 366	
	twoway ///	
		(scatter outcome_bergen_commuter2 ttime ///
			if xaxis >= $lastyearpre ///
			, connect(l) sort lcolor(black) mcolor(black) lpattern(dash) ///
			msize(small) yaxis(1 2)) ///			
		(scatter outcome_bergen_commuter ttime ///
			, connect(l) sort lcolor(black) mcolor(black) msize(small) /*lpattern(dash)*/) ///	
		, graphregion(color(white))	///
		xtitle("") ytitle("Electric vehicle share (pp)") ///
		name(figC2, replace) ///
		legend(on order(2 1) cols(1) ///
			label(1 "BEV share, paying commuters in Bergen") ///
			label(2 "Predicted BEV share in absence of the policy") ///
			pos(6) ///
		) ///
		xline(`announcement_date', lpattern(dot) lcolor(black)) ///mdy(2,1,2016) - implementation date
		xline(`implementation_date', lpattern(shortdash) lcolor(black)) /// mdy(2,18,2015) - announcement date
		ylabel(, angle(0)) ylabel(, angle(0) axis(2)) ///
		xlabel(2011(1)2017)
		
	graph export "${figures}figC2.png", as(png) replace
	graph export "${figures}figC2.pdf", as(pdf) replace
}
